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MODELING OF TURBULENT TRANSPORT 
IN THE SURFACE LAYER 

By G. Louis Smith 
Langley Research Center 

SUMMARY 

The turbulence equations as written by Donaldson using the method of invariant 
modeling have been applied to the following limiting cases of the surface or constant flux 
layer of the planetary boundary layer: 

(1) Neutrally stable 

(2) Stable (above influence of surface roughness) 

(3) Nearly neutrally stable 

(4) Very unstable (free convection) 

For the neutrally stable case, the equations are shown to admit as a solution the familiar 
logarithmic profile. By use of this result, boundary conditions suitable for the surface 
layer are defined and are simple to apply to rough surfaces. 

Expressions for the macroscale length A are given for each case. The parameters 
b, relating the microscale to the macroscale, and c, the ratio of macroscale length to 
height, are computed to fit atmospheric data. Owing to the structural difference between 
mechanically and thermally produced turbulence, different values are found for b and c 
for the neutrally stable and unstable cases. The b- value for the stable case agrees very 
closely with that for the neutrally stable case. 

INTRODUCTION 

The most studied case of atmospheric turbulence is that within the surface layer, 
which has a thickness on the order of 20 to 200 meters (ref. 1, p. 100). Many data have 
been taken by use of instrumented masts, and continue to be taken as more sophisticated 
instrumentation becomes available. These data have spurred the development of a body 
of theory, which may be used for analysis of data or for prediction of turbulent transport 
for given meteorological conditions. The basis of much of this theory is the Moran- 
Obukhov similarity theory (ref. 2). By this theory a characteristic velocity, length, and 
temperature can be defined such that all quantities when nondime ns ionali zed by using 
these three parameters are reduced to a set of universal curves, which have been empiri- 
cally defined (refs. 3 and 4). 



The surface layer is dependent upon the roughness of the surface beneath it, which, 
for example, may be grass, crops, trees, or buildings. This effect is accounted for by 
use of the roughness length, which is found empirically and is tabulated for various sur- 
faces (refs. 5, 6, and 7). 

One approach to the computation of atmospheric turbulence is that of Donaldson and 
others (refs. 8 and 9), who modeled second-order correlations of the Navier-Stokes equa- 
tions. Results are presented for application to the surface layer in reference 10. In 
reference 10, the boundary conditions suitable for a flat plate were used; that is, velocity 
and turbulence vanished at the surface. Also, the mean velocity and temperature profiles 
were specified (by using data from the work described in ref. 4), and the profiles of tur- 
bulent quantities were computed. Donaldson’s equations require the specification of a 
macroscale mixing length. This length is typically taken to be proportional to the distance 
from the surface out to some point, which must be selected, beyond which the macroscale 
length is assumed constant. 

This paper presents a study of the application of Donaldson s equations to the sur- 
face layer in light of existing surface-layer theory. This approach gives a guide to the 
application of Donaldson’s equations to atmospheric problems and conversely provides a 
possibility of generalizing surface-layer theory to a broader class of problems. The 
approach used in references 8 to 10 is to develop a single comprehensive theory of tur- 
bulence with a single set of constants, which will be applicable to all cases, from turbulent 
flow over a flat plate in a wind tunnel to atmospheric turbulence. The philosophy of the 
present paper is that atmospheric turbulence is sufficiently important to warrant compu- 
tation of a set of constants specifically for this case. In this paper, Donaldson’s equations 
are simplified for the large Reynolds numbers characteristic of atmospheric turbulence. 
Also, the turbulent Prandtl number is introduced into the model. The boundary conditions 
suitable for a rough surface (e.g., grass, forest, or city) are considered, and heuristic 
arguments are given for a set, which is selected. The equations are then solved for a 
neutrally stable atmosphere (no vertical heat flux), thereby resulting in the familiar loga- 
rithmic wind profile. This serves three purposes: The suitability of the governing equa- 
tions is verified; the suitability of the boundary conditions is verified; and two constants , 
in the method are evaluated for atmospheric application. Next the case of stable stratifi- 
cation above the influence of the surface roughness is treated. From this study comes a 
mixing length appropriate to the stable layer. This length can be used to determine the 
point above which the mixing length becomes constant. Finally, the free-convection or 
unstable case, in which turbulence is thermally produced, is considered. 
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SYMBOLS 


a constant relating macroscale of turbulence to microscale (see eq. ( 10 )) 

A n set of coefficients in series solution of temperature equation (see eq. (62a)) 
b constant relating macroscale of turbulence to microscale (see eq. ( 10 )) 

B constant in equation for free convection (see eq. (94)) 

B n set of coefficients in series solution of temperature variance equation (see eq. (62b)) 

c constant relating macroscale length A to height (see eq. (39)) 

C n set of coefficients in series solution of horizontal heat -flux density equation 
(see eq. (62c)) 

set of coefficients in free-convection solution (see eqs. (77) to ( 86 )) 

2 

g gravitational acceleration, m/ sec 

H constant for free-convection case defined by equation (100) 
k von Karman constant 

2 2 

K twice total turbulent kinetic energy per mass, + S 22 + S 33 , m /sec 
L Obukhov length (see eq. (23)), m 

L 2 length appearing in analysis of nearly neutral atmosphere (see eq. (63)), m 
n index in series solutions 

Pr turbulence Prandtl number 

q = T f (Vi - vi) (proportional to heat flux in i -direction), m-K/sec 
Q vertical heat flux (positive upward), kg-K/sec-m 2 

r temperature variance, T’ 2 , K 2 
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R a turbulence Reynolds number (see eq. (11)) 

s negative of turbulent shear stress, kg/sec 2 -m 

Sjj turbulent stress tensor, (vj - Vi)(Vj - Vj), m 2 /sec 2 

t time, sec 

T mean temperature, K 

T* temperature fluctuation from mean, K 

T Q mean temperature at z Q 

T* scaling temperature (see eq. (48)), K 

u mean velocity in direction of mean flow, m/sec 

u + friction velocity, m/sec 

Vj, velocity component in i-direction, m/sec 

w vertical component of velocity, m/sec 

x coordinate in direction of mean flow, m 

z coordinate in vertical direction, m 

z Q roughness height, m 

“ - Wfe/£ ) 

0 constant in solution form for free convection (see eqs. (72)) 

S = 1 + 2b -1 c 2 
3 

A macroscale length, m 
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A microscale length, m 

^ molecular viscosity, kg/ m-sec 

P 2 second viscosity coefficient, kg/m -sec 

v exponent in equation (72) 

o 

p density, kg/m 

y w > ctj variance of w and T, respectively 
^ nondimensional temperature gradient (see eq. (58)) 
nondimensional velocity gradient (see eq. (57)) 

The bar (") over a symbol indicates time averaged value. 


GOVERNING EQUATIONS 


A set of equations for turbulent flow applicable to atmospheric motions is given in 
reference 8, along with a suggested modeling of such correlation terms as is required to 
close the set of equations. In reference 9 these equations with this modeling are applied 
to the case of steady-state parallel turbulent shearing flow in the atmosphere over a 
large uniform plane. For this case, the flow is a function only of the vertical coordinate z 
and the continuity equation is trivially satisfied. With no body force or radiative heat 
source, the remaining equations reduce to 


Momentum: 


d u d 


dz dz 


(p S 13 ) = 0 


Energy: 


d 2 T d 


dz 


2 dz 


( P Qg) - 0 


( 1 ) 


( 2 ) 
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Longitudinal turbulent energy: 


-2 P S 


13 


du + A ( P VK A 


dz dz 


dz 


. P ^(s n .i K ) +li fh±-^s n = 0 

> 1 3 / d Z 2 x 2 


A 


( 3 ) 


Lateral turbulent energy: 


L fpvTCA^-p ^ fs 22 - i k) + p i ! 22 -i^s 22 = 0 

dz \ d z J A V 3 / dz 2 2 11 


(4) 


Vertical turbulent energy: 


dz 


(pvTcA^Vai^f. 

\ dz / dz V 


avk!!^ -pil IS 
dz / A 


33 


K 


-Pp-if 33 - liL S33 + 5fpq 3 - 2 (p + p 2 ) S 
dz 2 ,2 33 T 3 2 


33 


d 

(- 

dz 

\ p dz / 


£ (dp 

ri 2 \ dz 


= 0 


(5) 


Shear stress: 


-p S 


du 


dS 


■„ ~~ + 3iL p Vk a — 
^ dz dz V dz 


13 


- p S n + AVK^L 3 ^ 
A J dz dz 


+ M 


~ ~13 2p Q / A \ q 
— 7 - r s i3 “ ^ + ^ s i3 
dz ' 5 


d ( 1 dp\ 1 /dp 


dz dz/ p 2\dzy j 


pg„ 

[t 


+ qi = 0 


( 6 ) 
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Longitudinal heat flux: 


- pS 


13 


dT 

dz 


du 

pq 3 — + 

6 dz 


d U KA ^l\ . „ Pr VK „ d 2 qi . 2 m Q 
Iz \ Pr dz y A dz 2 2 


(7) 


Vertical heat flux: 


- pS 


33 


+ 3 d / 'W KA . p Pr YK 


dz 


dz \ Pr dz 


q 3 +P 


j2 

d q 3 


dz' 


2 p 


^3 


+ £§ir 



( 8 ) 


Temperature variance: 

-2Pq„5T + d (p^5E\ + "— -l5r = 0 
dz dz \ Pr dz/ dz 2 x 2 


(9) 


In these equations appear two lengths, A and X 3 which are the macroscale and microscale, 
respectively, of the problem. For turbulent transfer of heat terms in the present paper, 
the mixing length is divided by the turbulent Prandtl number Pr, The introduction of this 
modification to the original model of references 8 to 10 is in analogy with laminar trans- 
port equations. In the model of references 8 to 10, it is assumed that 


x = 


A 

V a VbR^ 


where 

R - ^ A 
A it 


( 10 ) 


( 11 ) 


Equation (10) was assumed in references 8 to 10 because it behaves for large and small 
r as desired and permits self-similar solutions of the decay of a free jet. Also, a = 2.5 
and b = 0.125 were selected in references 8 to 10, based on extensive wind-tunnel data. In 
the present paper, the value of b applicable to atmospheric problems is considered. For 
atmospheric processes, turbulence Reynolds number bR^ » a, and R^ is very high so 
that 

iiL = 2b gff?- (12) 

X 2 
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Where A appears as a mixing length for temperature or a correlation involving temper- 
ature, it is divided by Pr. Equation (12) is substituted into equations (3) to (9), and the 
remaining terms containing but not containing \ , are neglected as these terms are 
for molecular transport, which is much smaller than turbulent transport in the atmosphere. 
The momentum and energy equations reduce to 


— ^ S 13^ ~ 0 
dz 


-i (pq 3 ) = 0 
dz 15 


and are integrated to yield 

pSjs = Constant = -s = -pu* (13) 

P q 3 = Constant = Q (14) 

where u^. is the friction velocity a These are the familiar results that shear stress and 
heat flux are constant through the surface layer. Gradients of density are neglected. 
Equations (3) to (9) now reduce to 

Temperature variance: 

d / P ' / K A dr\ - 9V> Pr/>VK" r = 2 q51 (15) 

di ( ~pT dij A dz 


Heat flux: 


x: 


d (pVlCA dq A . (1 + 2b) ?rp>/g 


dz V Pr dz 


A 


-s^ + Q?!i 

dz dz 


(16) 


z: 


+ 2b)Q = pS^n^I - ££: 
A 66 dz T 


(17) 


Covariance tensor: 


'll- 


A (pVKA^H 

dz \ dz , 


- 2b S u = JL 
A u 


^ (S u -M - 2s^ 

A V 11 3 / dz 


(18) 
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( 19 ) 


22 ' 


i 


dz 


dz 


- 2b ^ S 22 = ^f( S 22 


K 


S 33 : 


i_ (pVK A ^ll) - 
dz \ dz / 


2b 


.^c, pVK 


33 


A 


S 33-J K 


_ 2gQ 


(20) 


*13’ 


pS„ *L - (1 + 2b) ^ g . Pg = o 
dd dz A T 


( 21 ) 


BOUNDARY CONDITIONS 

In classical aerodynamics, flow over a smooth surface is considered for which the 
boundary conditions are that all components of velocity vanish. This applies not only to 
the mean velocity, but as well to the turbulent part. The vanishing of the turbulence at the 
surface requires that momentum transport be carried by a laminar sublayer. For the case 
of the atmospheric surface layer, the surface may be smooth, as a snow-covered plain, 
but is typically rough, as a forest, a city, or crops in a field. Although each leaf in the 
forest may have a laminar sublayer, there is a lack of enthusiasm for considering this 
level of detail in specifying the geometry. Also, much of the momentum transfer at the 
rough surface is caused by the drag of the wind on trees or houses. 

For a neutrally stable surface layer, that is when Q = 0, many studies have shown 
that the surface-layer wind profile above the roughness elements may be described by a 
logarithmic function 



( 22 ) 


where z 0 is the roughness height characterizing the surface (refs. 1, 5, and 6). This 
relation is written such that u = 0 at z = z 0 . The roughness height z Q is tabulated for 
different types of surfaces on page 233 of reference 5 and page 150 of reference 6. 

The boundary conditions for the turbulence terms S^, S 22 j S 33 , and r are arrived 
at by use of the Monin-Obukhov similarity theory together with data. In references 1 and 3, 
it is pointed out that according to Monin-Obukhov theory, S^/Sjg, S^/S^, S 33 /S 13 , and 
r/T* are functions only of z/L, where L is the Obukhov length given by 


3 

u£ p T 

L = - 

kgQ 


( 23 ) 
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and T + is the scaling temperature given by 



Because u + and Q are constant through the surface layer, L and T* are each constant. 
These functions are defined from data in chapter 4 of reference 1 and in reference 4, and 
it is pointed out that S.^, S 22 > and vary very slowly with altitude and are moreover 
constant for the neutrally stable surface layer. On this basis, the boundary conditions, 
which are proposed in this paper for these quantities at the rough surface, are 


dS 


11 


dS 


22 


dS 


33 


dz 


dz 


dz 


= 0 (z = z 0 ) 


(24) 


These conditions state that near a rough surface turbulence is generated at the same rate 
as it is dissipated. Use of these conditions assures that the surface layer will approach 
neutrally stable layer behavior. 

From the discussion of temperature variance in chapter 4 of reference 1, \/r/T* 
is constant with altitude for an unstable surface layer for small z/L and for a nearly 
neutral surface layer. Also, for small z/L in a stable surface layer, d(-/r/T*)/dz = 0. 
Thus the boundary condition to be used herein is 

~ = 0 (Z = z 0 ) (25) 

dz 

Equations (15), (16), (18), (19), and (20) are second order and as such require two 
boundary conditions each. The form of these equations, taken one at a time, is such that 
there is a homogeneous solution which vanishes exponentially with altitude and one which 
grows exponentially with altitude. For these equations, the exponentially increasing term 
is eliminated. A situation where turbulence is generated above the surface layer can be 
conceived. This turbulence would be damped exponentially as it diffused downward. By 
using boundary conditions, which get rid of the exponentially increasing part of the solution, 
this situation is omitted. 

One more boundary condition needed for the longitudinal heat-flux equation is that 
the q^ -profile be well behaved at the rough surface. 

Note that, had viscous terms been retained, the surface would not be a singular point 
of the differential equations governing q^. There would then be no occasion to require that 


10 



the solution be well behaved, but there would be the laminar boundary condition that 
= 0 at z — zq . 


NEUTRALLY STABLE ATMOSPHERE 

The neutrally stable atmosphere is now considered. The purpose of this is not to 
present another derivation of the familiar logarithm profile but to verify the application 
of the equations of reference 8 and the use of the boundary conditions proposed here for 
analysis of the surface layer. Also, a reevaluation of some arbitrary constants in the 
model is desired for atmospheric use. 

The neutrally stable atmosphere is defined by the vertical heat flux Q = 0. A solu- 
tion to equations (15) to (17) is that r, q, and dT/dz vanish throughout the surface layer. 
Equations (18) to (21) are solved by assuming that throughout the surface layer 


dS ll dS 22 dS 33 


= 0 


dz dz dz 


so that these equations reduce to 


■ZbpVK S ni = pVK (S n1 - - K| - 2s A ~ 
11 1 11 3 / dz 


■2b p s/K S 22 = K (S 22 - ± K 


-2b pVK S oo ” /? VlC [S o o *" — K 


>33 


'33 


Summing equations (27) to (29) and substituting K ~ + S22 + S33 yields 


b*K 3/2 = SA* 
dz 


Equations (28) and (29) give 


S„„ = S 


K 


22 - 33 “ 3(1 + 2b) 


(26) 


(27) 

(28) 

(29) 


(30) 


(31) 
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From equations (27) and (30) 


11 


1 + 6b 
3(1 + 2b) 


K 


Equations (30), (31), and (21) are combined to produce the result 


K = ±{1 + 2b) — 

b p 


( 32 ) 


(33) 


The friction velocity u* for turbulent flow is defined by 


2 

pu* = s (34) 

Equations (31), (32), and (33) thus yield 


S 22 S 33 1 

<4 4 " ^ 

(35) 

S 11 1 + 6b 

«T ^ 

(36) 

K 3(1 + 2b) 

u 2 VSE 

(37) 

Equation (30) becomes 


A — = 3 3 ^ 4 b 1//4 (l + 2b) 3 / 2 u* 
dz 

(38) 

At this point, the mixing length A must be specified, 
form assumed for A is 

In reference 8 and elsewhere, the 

A = cz 

(39) 


Thus the solution to equation (38) is. 
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u = 


r 3 3 / 4 b X / 4 (1 + 2b) 3 / 2 



Note that the roughness height is thus incorporated into the solution, as typically is done 
in surface -layer theory. Comparison of this solution with equation (22) shows that 


3 3/4 b l/4 (1 + 2b) 3 / 2 = 1 ^ (40) 

c k 

The results obtained are consistent with the initial assumption of equation (26) and satisfy 
the governing equations and boundary conditions; hence, the solution is self consistent. 

The theory is now compared with empirical results. The figures in reference 3 show 
that for a neutral atmosphere 



2.6 



=: 2 


From these values 

*L = 12.45 
u* 

Equation (37) may now be solved for b, thus giving b = 0.021. Equations (35) and (36) next 
give 




2.12 
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= 2 


Vs 


33 _ 


^ 2 


u* 


By assuming von Karman's constant k = 0.35 for atmospheric application (ref. 4), equa- 
tion (4) gives 


c = 0.323 

As applied, the present model is seen to predict a longitudinal turbulent energy which is 
lower and a vertical turbulent energy which is higher than data. Fortuitously, the lateral 
turbulent energy is in agreement. The arbitrariness of the results is apparent; equa- 
tion (35) or (36) could have been used to evaluate b, thereby changing the numerical 
results- Equation (37) is solved for b because the total -turbulent-energy term K is felt 
to be more representative than any one of the three parts Sjj, 822? or S33. Also, K is 
the only one of these quantities to appear in each equation. 

Thus far, it has been shown that the model developed by Donaldson, Sullivan, and 
Rosenbaum where used with the boundary conditions previously discussed does describe 
the neutrally stable boundary layer, if the partitioning of the turbulent energy among the 
longitudinal and vertical components is of no concern. This model uses a single micro- 
scale length X to describe the dissipation of turbulent energy. In reference 11, Donaldson 
begins with three microscale lengths to describe the dissipation. In reference 8, the 
microscale length relates one second-order tensor to another. The simplest choice is to 
let the microscale length be a scalar, which permits only a single length. In order to 
incorporate more than a single length in the modeling of the dissipation term of the velocity 
correlation, the microscale lengths must be modeled as components of a fourth- order 
tensor. To do so is lengthy and is not done here. 

The values derived herein for b and c differ considerably from the values b = 0.125 
and c - 1.58 found in reference 8. The reason for this disagreement is that the present 
paper uses atmospheric data to evaluate these numbers and reference 8 uses flat-plate 
data from reference 12. 
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STABLE ATMOSPHERE 


A stably stratified surface layer, that is, one in which the temperature increases 
with height, is described quite well by a log-linear profile. Near the ground, the logarith- 
mic term dominates. Well above the ground the velocity and temperature gradients are 
nearly constant, as are the turbulence and mixing length A. For this condition, the 
governing equations (15) to (21) simplify considerably, thereby reducing to 

Temperature variance: 

-2b = 2Q~ (41) 

A dz 


Heat flux: 


x: 


z: 


Covariances: 


S 


S 


S 


S 


11 : 

22 ' 

33* 

13 : 


X 9V0 dT ~ du 

(1 + 2b) — q-, = - s + Q 

A dz dz 


-(1 + 2b) Q Pr '{K = PS 33 - £l 


dz T 


-2b^S 11 =^I(s 11 -lK) -2s^ 

A 11 A V 11 3 j dz 


-2b^5s 22 =^ (s 22 -iK 


-2b 


P^K « __ P / 1 

~X- b 33 "~A~ l b 33 “~ 


K ) - 
3 / T 


pSoo ~ - (1 + 2b) s - £* q- = 0 
66 dz A T 1 


(42) 

(43) 

(44) 

(45) 

(46) 

(47) 


Equation (44) shows that Sjj is driven by the wind shear du/dz, and equation (46) shows 
that S33 is reduced by the downward heat flux. In this section, quantities are nondimen- 
sionalized by use of Obukhov length, scaling temperature, and friction velocity. Nondi- 
mensional quantities are denoted by a caret ('). Thus 
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du _ L du 
dz u * dz 



U* 


dT = L_ dT 
dz T* dz 



u^T.^ 


J 


(48) 


According to similarity theory, the nondimensionalized quantities may be expressed as 
functions of nondimens ional height z = z/L. 

Equations (41) to (47) are nondimensionalized in this manner and become 
Temperature variance: 


r = k 


dT A 
d z b Pr fk 


Heat flux: 


x: 


(1 + q, Pr = 

A 1 dz dz 


z: (1 + 2b) kIKpr = S„ £1 - JL 

A dz k 2 


Covariances: 


c . A _ K 2 A du 

ir 11 ~WTzii> £ a ; 2b) di 


S 22 : S 22 


K 


3 (1 + 2b) 


(49) 


(50) 

(51) 


(52) 


(53) 
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2 A 


(54) 


S 33 : S 33 = 


K 


3(1 + 2b) p + 2b) kjk 


S 13 : S 33 


du 

dz 


- (1 + 2b) 


KA = 0 


A k 2 


Turbulent kinetic energy: 



A 


du 

dz 


_1 

k 


(55) 


(56) 


Equation (56) is obtained by summing equations (52), (53), and (54). 

Here are eight equations in 11 unknowns: K, S^, §225 S 33 , du/dz, dT/dz, q^, A , f, 
Pr , and b. Three additional requirements may be imposed on the problem. Two require- 

^ ^ -A 

ments are that du/dz and dT/dz match empirical values. A third requirement is that K 
match the empirical value. 

As before, K is derived from the figures in reference 3; thus 


n/S 


11 _ 


= 2.3 


u, 


or 


and 


V ^22 


= 2 


U* 


vS 33 

= 1.3 

u* 


?11 = 

5.29 

^22 = 

4 

^33 = 

1.69 


K = 10.98 

Experimental studies (ref. 4) have shown that nondimensional wind shear 

t-kz*" 
m d 2 
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and nondimensional temperature gradient <t> = z _ are related to z by 

n a ™ 


<}> = 1 + 4.7z (S'?) 

m 

0 h = 0.74 + 4.7z (58) 

For z » 1, that is, above the influence of ground roughness, these equations result in 
empirical values of 

k^=4.7 (59) 

dz 

£1=4.7 (60) 

dz 

For K = 10.98 and k = 0.35, equations (59) and (56) give 

JL = 0.29 

A 

thus equation (49) gives 


Pr r = 1.71 


From figure 4.23 of reference 1, r « 0.09 for the stable layer, so that for Pr * 0(1), the 
present equations result in a value of r, which is high by a factor of 20. The solution to 
the set of equations (50) to (55) using these results is found to consist of complex numbers. 
Equation (49) was replaced by r = 0.09, and equations (50) to (55) were solved, thereby 
giving 


S 


11 


4.12 


§22 - 3.50 


S 33 = 3.37 
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b = 0.0232 


A = 0.0798 


Pr = 0.99 


= 0.2185 


The value of b thus computed for the stable case is quite close to the value of b computed 
for the neutrally stable case. This agreement increases confidence in the method. The 
comparison of the computed values of S 11# § 22 , and S 33 with the measured values shows 
the agreement to be fair for and § 22 , but rather poor for Because part of the 
S 22 measured is due to small-scale horizontal motions as well as mechanical turbulence 
and heat convection (ref. 3), the computed § 2 2 should be lower than the measured value, 
as it is.. 


Although the differences between the b-values computed here for the neutrally 
stable and the stable atmospheres are minute, a single b-value is desirable. This value 
was arrived at by slightly reposing the problem as a least squares fit. For a given b, the 
K, S u , S 22 , and S 33 were computed for the neutral and for the stable atmosphere, and 
the sum of the squares of each of these quantities from its measured value was computed. 
The b for which the sum of the squares was minimum was b = 0.0222. This value is 
considered as the best fit for the two cases jointly. The corresponding mixing length in 
the stable case is A = 0.0785. 


In reference 4, a - 


9 3 du/dz 


is found to be approximately 1.2. Comparison of 


S i3 dT/dz 

the definition of a with that of Pr shows that Pr = 1/a so that empirically Pr = 0.83 for 
the stable case. Thus, good agreement of the vertical heat-flux equation with data is 
shown. The factor of 20 between the measured and the calculated values of r indicates 
that additional work is needed in the modeling of the temperature -variance equation. 


NEARLY NEUTRAL ATMOSPHERE 


The case is now considered for which the heat flux Q is small and the influence of 
Q on the turbulent motions is small. Inspection of equations (18) and (20) shows that this 

case requires that 

glQL< (61a) 

T dz 
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and from equation (21) 


Pg 

T 


N « 



(61b) 


For these conditions, S„ is essentially unaffected by Q and is given by the solution for 
the neutrally stable ground layer and A is given by equation (39). Equations (15), (16), 
and (17) can now be solved for r, q^, and dT/dz by assuming series solutions 



Here 

PS 33 TVK 
2 g C Q Pr 


(62a) 


(62b) 


(62c) 


(63) 


is a parameter with the dimension of length, which appears in the analysis, and is quite 
similar to the Mo nin- Obukhov length. Equations (62a) and (62b) are substituted into 
equation (15) so that the are expressed in terms of the A^. Thus 


2S 33 \ A n 
, c 2 K ) n 2 c 2 - 2bPr 2 


(n = 0, 1, 2, . . .) 


(64) 


This result is used in equation (17) so that the are found to be given by 
Aq = -(I + 2b) 


(65a) 
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n+1 


n 2 c 2 -2bPr2 


(n =0, 1, 2, . . .) 


(65b) 


The temperature gradient is thus 

dT = _ Pr V K Q ( 1 + 2b) 
dz 


PS33 c 


z (1 + 2b)L 2 ^ 


CO 

L *n 


z 


n-1 


( 66 ) 


and the temperature profile by integration of equation (66) is 


T = T 0 - 


Pr'/KQ (1 + 2b) 


pS 33 c 


"in z - — _ — E ^(4 

(1 + 2b) n=l n l L 2 J 


(67) 


Equations (64), (65), and (62b) give 


= Lll b + 2Pr 2 f 

2 q \ b “j 

p Sgg \ n-1 


A (JL 

n [ T 

z 


1 2 2 „2 

n=l n c - 2bPr 


( 68 ) 


Equation (16) can now be solved for the horizontal heat flux, thus giving 


u* Q 


*1 


PS 


33 


-1 - 


cS 


33 


k(l + 2b) Pru*v1C 
A. 


+ Pr 


n=l e 2 n 2 - (1 4- 2b)Pr 2 \ L 2, 


00 2 

For nc »2bPr (orn » 1), equation (65b) results in 


A a 

n o 

(n!) 2 


(69) 


so that the series in equations (66) to (69) converge absolutely for all z, and quite rapidly 
for computations. However, the solution was based on the approximation that heat flux 
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does not influence the turbulence, as expressed by inequalities (61a) and (61b). By 
equation (69), jq-jJ % u*| Qj /,> vK; therefore, with the use of equations (37) and (38), it 
is seen that the inequality (61a) implies (61b); thus 



s|q| b|Q| 

T VK * T 



s]q! 

T 


« s 


du 

dz 


f'j 


s 


Vk 

A 


Equations (30), (35), and (37) are used to rewrite the right-hand side of inequality (61a); 
thus 


g 1 Q ( _ 3b > ls/K S 33 
T cz 

This result is rearranged and equation (63) used, thereby leading to the condition 
for validity of the solutions of this section that 

— « 3b 

N 


For this reason, only the leading terms in equations (66) to (69) need be retained. 
From equation (66), for z/Lg — * 0 


dT 

T* dz 


r nr 

V K k (1 + 2b) 


S 33 c 


Pr = Pr 


(70) 


where equations (35), (37), and (40) have been used to reduce the expression within the 

brackets to unity. Measurements of the surface layer (ref. 4) show that for z/L 9 — * 0, 

z dT " 

— * 0.74. Thus, Pr = 0.74 for the nearly neutral surface layer. This result 

T ^ dz 

agrees fairly well with the value of 0.99 derived for the stable case. From equation (68), 
for z/Lg — ■> 0 


Q 
r = 


1 + 2b 


P 2 s 


33 




’33 


1 + 2b 


(71) 


o 

From reference 1, measurements show that r/Tj * 1 for a nearly neutral atmosphere. 

For b = 0.0222, equations (71) and (35) give a value of t/t| - 1.49, which is fair agreement. 
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FREE CONVECTION 


The case is now considered of a surface layer with a strong upward flux of heat 
and low shear. In this condition, called free convection, the turbulence is primarily 
thermally produced (ref, 1) in contrast to the condition in which turbulence is produced 
by shear. Thus, in equation (18) the term s ^ is neglected, and the turbulence is 

Q2 

driven by the term gQ/T in equation (20). 

A solution form 

VK = 0z v (72a 

A = cz 1 ' (721 

is assumed, and a solution to equations (18), (19), and (20) is found to be 

A = cz (73 


VK = F, (_sQ\ 1 ' /3 z l/3 


"a ■ A ”2 (f r « 2/S 


S 11 = S 22 


3 ( 1 + 2b - - c 2 


gQ \ 2/3 z 2/3 


where 


1 + 2b - — c 2 ) F ? - i 
3/^3 


F 2 = l- 


3^1 + 2b --c 2 


With A and K determined, equations (15) and (17) may be solved simultaneously for 
dT/dz and r, thereby resulting in 


= - F 3 ~ 


T\ 1/3 /Q\ 2/3 ,-V3 


f Z 
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* - r.gf er*- !/3 

where 

(80) 

EE (1 + 2b)Fj 
„ c 

f 3 

3 c 

F 1 F 2'^ 

9Pr - bPr 

(81) 

Fo 

f 4 = C — 

C 2 Fj 

bPr - 1 

9Pr 

(82) 

The shear du/dz and the horizontal heat flux qj 
tions (16) and (21). The result is 

can now be computed by solving equa- 

*i. f 

dz & p \Qg / 

(83) 

- ■ *.ensr-** 

(84) 

where 


Fc = . 1 + 2b + -A- 

c FjFj f 2 

(85) 

1 + . 2b + F, 
cF 1 F 2 

(86) 

6 „ nj c (i + 2b) Pv _ i - 1 

‘•s' . F?Fo 

- "_J 


This solution is based on the assumption that the turbulence generated by buoyancy is 
much larger than that from shear, so that 

du g M 

s — « B.Q 

dz T 

(87) 
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Equation (83) may be used to get the result that 


F 5(T 


« 


gQz\ 

T P / 


4/3 


for the free-convection solution to be valid. By use of the Obukhov length L, equation (23), 
this condition may be written as 



« 



( 88 ) 


The solutions of this section do not show the arbitrary constants of integration. In obtain- 
ing the self- similar solution, it is intrinsically assumed that the effect of the lower 
boundary has damped out at the altitude required to satisfy condition (88), and that like- 
wise the upper boundary is sufficiently remote that its influence does not diffuse into the 
region of interest. Thus the terms containing these constants are neglected. 

The solutions of this section may be determined by dimensional analysis except for 
the constants (refs. 2, 13, and ch. 3 of ref. 1). For this case, u* is negligible, and the 
governing variables reduce to Q/<o, z, and g/T. 

If the values b - 0.0218 and c = 0.323 computed for the neutrally stable case are 
used in equation (74), vie is computed to be negative. The possibility that \/K may be 
negative is eliminated by consideration of equation (20) for S 33- It is concluded that b 
and c must be quite different for thermally produced turbulence from their values for 
mechanically produced turbulence. When the structure of the turbulence for the two cases 
is considered, this difference is not surprising. For the mechanically produced turbulence 
of neutrally stable or stable turbulent flow, the turbulence appears as eddies, which have 
essentially a two-dimensional rolling motion in the direction of the flow. For free con- 
vection, the turbulence elements have the familiar mush room -cloud configuration, being 
axisymmetric about a vertical axis. Mechanically produced turbulence entrains fluid 
along a sheet, which may be visualized as being rolled up by the eddy. For thermally 
produced turbulence, fluid is entrained along a ring at the base of the mushroom- cloud 
top. These concepts are shown in the following sketches: 
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Fluid 

entrained 



Mechanically produced 
turbulence 



Fluid 

entrained 


Thermally produced 
turbulence 


Because of this distinction in structure, thermally produced turbulence is very effective 
in transferring heat, but not momentum. Reynolds analogy does not apply here (ref. 14); 
therefore, the turbulent Prandtl number may be quite different. The computation of a 
suitable b, c, and Pr is now considered. 

There are some basic requirements for b and c. First, the mixing length A must 
be positive; thus 

c > 0 (89) 


Next, = S 22 > which by use of equation (76) implies 


b 


1 2 

> — c 


1 

2 


(90) 


A more stringent constraint is given by the requirement that > 0, which by use of 
equations (75) and (78) implies 


b 



1 

6 


(91) 


Finally, it has been noted that VK > 0. This requirement implies, by equation (74), that 
F^ > 0; thus after some manipulations 


c 2 < - 5 " 

8 S 

3 

where 


5= 1 + 2b -lc 2 
3 


(92) 


(93) 
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These constraints are shown in figure 1. The condition (92) is applied by determining 
the boundary in the 5 -c plane, and then by mapping this line into the b-c plane of fig- 
ure 1, which shows the regions of permissible values of b and c. It is seen that for c>0, 
the requirement that Vlf > 0 is more stringent than the requirements that S^, S22? 
and Sgg be positive. Figure 1 also shows that the point b,c for mechanically produced 
turbulence lies well outside the region of permissible values for purely thermally induced 
turbulence. 

Three numbers, b, c, and Pr are to be computed from experimental data. Many 
possible conditions may be used. The conditions used here are that the theory and exper- 
imental data match for (1) vertical component of turbulent energy, (2) temperature 
gradient, and (3) temperature variance. These three criteria were chosen because they 
have been most thoroughly investigated (e.g., ref. 15). Alternate criteria would be to 
match the turbulent kinetic energy K or the velocity gradient du/ dz . 

The vertical component of turbulent energy S^g is given by equation (75), which 
may be rewritten as 

S33 = B 2 (IS ) 273 z 2 / 3 (94) 


where 


F 


2 

1 



(95) 


This notation agrees with that of reference 1. (Note that S^g -o w 2-.) Also, this form 
agrees with that of reference 15, with the notational change a2 = B. From reference 15, 
B = 1.4 (using k = 0.35) is found to be the value best fitting experimental data. With this 
number for B, equation (95) constitutes an implicit relation between b and c, which is 
shown as a curve in figure 1. 

The matching of the temperature gradient is now considered. In reference 4, the 
nondimensional temperature gradient <p ^ is shown for very unstable conditions. This 
gradient can be written approximately as 


\ - 0.23 



(96) 


Comparison of equation (96) with equation (79) shows that agreement between these two 
equations requires that 

F 3 = 0.23k" 4 / 3 = 0.93 (97) 
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Because Fg is a function of b, c, and Pr as given by equation (81), equation (97) is an 
implicit relation among these three parameters. 

Finally, the temperature variance is considered. Comparison of equation (80) with 
the experimental correlations of reference 15 for temperature variance shows that 

F 4 = a| , (98) 

By using values from reference 15, the right-hand side of equation (98) is found to be 1,63. 
As before, equation (98) provides an implicit relation among b, c, and Pr. 

Equations (95), (97), and (98) may be solved for b, c, and Pr, thereby resulting in 


b = 0.39 
c = 0.52 
Pr = 0.56 


The values for b and c are plotted in figure 1. It is seen that the values for b and c 
differ drastically from the neutrally stable cases. Also, Pr is considerably less for the 
unstable case than for the stable case, corresponding to the fact that the unstable atmos- 
phere is far more efficient at transferring heat than momentum. These differences from 
the neutrally stable and stable cases are attributed to the different structures of mechani- 
cally produced and thermally produced turbulence. 

As a check on this set of parameters, values from other sources were used. Chap- 
ter 4 of reference 1 quotes B as being tentatively 1.7 and shows data indicating that 


JL dT - 1 (99) 

° T dz 

for the unstable regime. In reference 13, Priestley gives the result 

H st 9 = 0.81 (100) 



Thus, from equation (100) it follows that 

F 3 = H- 2 / 3 (101) 
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which is in the form of equation (97). Equation (99) leads to the relation 

F 3 2 _ / z dT\ 2 = x 
F 4 \ ° T dz / 


(102) 


Numerical solution of equation (95) with B - 1.7 together with equations (101) and (102) 
leads to the result 

b = 0.45 
c = 0,59 
Pr = 0.65 

The values for b and c are plotted in figure 1. 

In concluding this discussion of the unstable case, it is pointed out that the applica- 
bility of the self-similar solution - — oc z“ 4 / 3 has been questioned (ref. 4). In refer- 

dz 

ence 15 the vertical velocity and temperature variances are found to agree with the free- 
convection form. In order for the self- similar solution to be applicable, the upper bound- 
ary must be high enough to insure that the motion is not influenced. 

CONCLUDING REMARKS 


The turbulence equations as written by Donaldson using the method of invariant 
modeling are applied to the following limiting cases of the surface or constant flux layer 
of the planetary boundary layer: 

(1) Neutrally stable 

(2) Stable (above influence of surface roughness) 

(3) Nearly neutrally stable 

(4) Very unstable (free convection) 

For the neutrally stable case, the equations are shown to admit as a solution the familiar 
logarithmic profile. By use of this result, boundary conditions suitable for the surface 
layer are defined and are simple to apply to rough surfaces. 

Expressions for the macroscale length A are given for each case. The parameters 
b, relating the microscale to the macroscale, and c, the ratio of macroscale length to 
height, are computed to fit atmospheric data. Owing to the structural difference between 
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mechanically and thermally produced turbulence, different values are found for both b and 
c for the neutrally stable and unstable cases. The b -value for the stable case agrees 
very closely with that for the neutrally stable case. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va,, June 11, 1973. 
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